Short-time Spin Dynamics in Strongly Correlated Few-fermion Systems 
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The non-equilibrium spin dynamics of a one-dimensional system of repulsively interacting fermions 
is studied by means of density- matrix renormalization-group simulations. We focus on the short-time 
decay of the oscillation amplitudes of the centers of mass of spin-up and spin-down fermions. Due 
to many-body effects, the decay is found to evolve from quadratic to linear in time, and eventually 
back to quadratic as the strength of the interaction increases. The characteristic rate of the decay 
increases linearly with the strength of repulsion in the weak-coupling regime, while it is inversely 
proportional to it in the strong-coupling regime. Our predictions can be tested in experiments on 
tunable ultra-cold few-fermion systems in one-dimensional traps. 



I. INTRODUCTION 

In an electron liquid the motion of one of the two spin 
species, e.g. in the presence of a spin current, can drag 
along the other one because of electron-electron interac- 
tions. This is the spin Coulomb drag effect or simply the 
Spin Drag (SD) [1-3]. In electron transport SD can be 
described by a frictional force proportional to the differ- 
ence between the velocities of the two populations and 
is described by a damping term in the equation of mo- 
tion for the time derivative of the spin-resolved center- 
of-mass momentum. SD has been observed [4, 5] in two- 
dimensional electron gases in semiconductor hetero junc- 
tions. 

The concept of SD can be extended to other quan- 
tum fluids with distinguishable species that can exchange 
momentum due to mutual collisions. Ultracold atomic 
gases [6] are clean systems in which SD can be observed 
in a truly intrinsic regime [7-11]. Further, the inter- 
action strength between atoms can be tuned at will by 
employing Feshbach resonances [6]. 

This work is motivated by a recent pioneering exper- 
iment [12] on SD in an equal mixture of two hyperfine 
states of 6 Li atoms confined in a trap. The authors 
of Ref. 12 measured independently the time-dependent 
position of the centers of mass of "spin-up" and "spin- 
down" particles starting from an initial condition in 
which the two types of particles are grouped in well- 
separated clouds. The experiment is performed in the 
"unitarity limit" in which the strength of interactions is 
the largest possible. At long times the separation of the 
centers of mass decays exponentially to zero. By measur- 
ing the time constant of this exponential decay the SD 
coefficient is determined. 

Besides providing information on SD in the strong cou- 
pling regime, Ref. 12 provides a wealth of new data on 
the short-time behavior - long before the SD regime is 
attained. There it is found that the two clouds perform 
several cycles of oscillation before settling at the bot- 
tom of the trap. If interactions are sufficiently strong, 
they reflect off each other several times before the inter- 



diffusion process begins. This short-time regime of spin 
dynamics, the short-time SD (STSD), constitutes the fo- 
cus of the present work. We tackle it non- perturb at ively 
by the time-dependent density-matrix renormalization 
group (TDMRG) method [13] (see Sec. A). This method 
is essentially exact, its main limitation being the maxi- 
mum system size that we can handle [14]. Starting from 
an initial condition similar to that of Ref. 12, we find 
that the oscillation amplitudes of the centers of the spin 
clouds decay in time quadratically for weak interactions, 
linearly for intermediate interactions, and again quadrati- 
cally for very large interactions. Below we argue that this 
intriguing reentrant behavior is a many-body effect. Our 
predictions are amenable to experimental testing, since 
in a recent work Serwane et al. [15] were able to trap 
few fermions in a ID geometry and to tune their mutual 
interactions by means of a Feshbach resonance. 

II. THE MODEL 

We consider a two-component fermion system with re- 
pulsive short-range interactions in a ID trap. The system 
is prepared in the ground state of two (spin-dependent) 
displaced harmonic potentials [Fig. la)]. At time t = + 
these external potentials are suddenly turned off and the 
system evolves in presence of a single harmonic confine- 
ment, according to the Fermi-Hubbard (FH) Hamilto- 
nian, 

H = - J^(ct a Ci + i, ff +H^O + ^S^t^a+S Wifli • 

i,a i i 

(i) 

Here J is the inter-site hopping parameter, c\ a (c^ a ) 
creates (destroys) a fermion in the i-th site (i G [1,£], L 
being the total number of lattice sites), a =t, I is a label 
for a pseudospin-1/2 (hyperfine-state) degree of freedom, 
U > is the on-site repulsion, = c\ a Ci y(T is the local 
spin-resolved number operator, and hi = n^ + n^. The 
third term on the r.h.s. of Eq. (1) represents an external 
parabolic potential Wi = Vi{% — L/2) 2 of strength V2, 
corresponding to a frequency uo = 2vT/^J /h. 
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We follow the time-evolution of the spin-resolved den- 
sities (hi j(7 (t)) on a time scale much smaller than the 
spin equilibration time [16] and calculate the spin- 
resolved centers of mass from Xcm a(t) = L~ x J2 i=1 (i — 

L/2)(^(t)\h^mt)). 

III. NUMERICAL RESULTS AND DISCUSSION 

In Fig. 1 we illustrate the time evolution of the oc- 
cupation numbers rii y(7 (t) for a system of N = 6 spin- up 
particles and N = 6 spin-down particles, in a lattice with 
L = 240 sites. The harmonic potential has a strength 
V2/J = (1/160) 2 , corresponding to a harmonic oscillator 
length aho = (J/V2) 1 / 4 ~ 12.65, in units of the lattice 
constant, and to a frequency uj = J/ (80ft). These pa- 
rameters have been used also for all other plots and their 
choice yields minimal lattice effects (see below) [17]. The 
data in Fig. 1 correspond to U/J — 5. In panel a) we il- 
lustrate the initial state, with two non-overlapping clouds 
with opposite spins. Panels b)-f) show the time evolution 
of this initial state. We highlight two features: i) in pan- 
els b) and e) high-density regions form near the center of 
the trap due to strong repulsive interactions [12]; ii) in 
panels c), d), and f) we see how the spin- up cloud (blue 
curve) drags along a substantial fraction of down-spin 
atoms (red curve). 

In Fig. 2a) we show the time evolution of spin-resolved 
cent er-of- mass, XcM,a(t), in the weak coupling regime, 
U / J < 0.05. In absence of the lattice, the center-of-mass 
of each atomic cloud is decoupled from "internal" degrees 
of freedom and should oscillate at the trap frequency, co>, 
without decaying. This is confirmed by the data corre- 
sponding to U / J = in Fig. 2a) (dotted lines). No visible 
damping effects appear within the time-scale of the plot, 
since we have minimized lattice effects [18]. 

When U/J is finite the two clouds still go through each 
other, but their motion is damped. Fig. 2b) reports the 
maxima of the blue and red curves as a function of time, 
for several different values of U/J < 0.05. The amplitude 
of the oscillations in Fig. 2a) decays quadrat ically in time. 
This is because, in this regime, the center of mass of each 
cloud is a harmonic oscillator weakly coupled to internal 
degrees of freedom. The relevant excitation spectrum, 
5(0), is sharply peaked about ft ~ uj. The position 
of the peak determines the frequency of the oscillations 
and the second moment of the spectrum determines the 
quadratic decay of their amplitude. The quadratic de- 
cay can be also verified analytically by means of time- 
dependent perturbation theory - see Sec. B. 

We now discuss the strongly correlated regime, U / J ^> 
1. The main results are summarized in Fig. 3. Dot- 
ted lines in Fig. 3a) represent the exact time evolution 
of the spin-resolved center-of-mass for U/ J = 00. In 
this limit, the Hamiltonian in Eq. (1) maps onto [19] 
Hoc = ~JVJ2i,a(4,a^i+i,o- + H.c.)P , where V is a 
Gutzwiller projector (that avoids double occupation of a 
lattice site). Dotted lines have been obtained by applying 
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FIG. 1: (Color online) Time evolution of the site occupa- 
tions ni,cj(t) for a system of twelve fermions at strong coupling 
(U/J = 5). Panel a) Initial state: two clouds of atoms with 
opposite spin (blue and red curves) are spatially separated 
using two displaced harmonic potentials. Panels b)-f) Subse- 
quent time evolution after the abrupt switch-off of these local 
potentials. The two clouds are forced to propagate against 
each other in presence of an overall harmonic confinement 
of frequency uj. The parameter t denotes time in units of 
the period T = 2tt/uj induced by the harmonic confinement. 
Dashed lines indicate the initial spin-dependent displaced har- 
monic traps [panel a)] and the overall harmonic trap for t > 
[panels b)-f)], and are plotted as guides to the eye. 



TDMRG to Hoc- Notice that the centers of mass of the 
two clouds behave practically like two classical particles 
that bounce off each other quasi-elastically oscillating at 
twice the trap frequency. The frequency doubling with 
respect to Fig. 2a) is understandable as follows. Due to 
strong repulsion, fermions of opposite spin are confined to 
one half of the trap: effectively, only the antisymmetric 
levels of the harmonic oscillator, whose energy separa- 
tion is 2cj, are involved in the time evolution. A rather 
complicated dynamical pattern, however, is present in 
the time evolution of the spin-resolved site occupations 
ni,<r(t) (see Fig. 1). 

In Fig. 3b) we plot the amplitude of oscillations vs 
time for 5 < U/J < 30. Remarkably they decay lin- 
early. As mentioned in the Introduction, the quadratic- 
to-linear crossover is a many-particle effect. One can in- 
deed solve analytically the evolution dynamics for an in- 
teracting system of two particles with antiparallel spin in 
a harmonic potential. In that case, the time evolution of 
XcM,a follows the quadratic behavior seen in Fig. 2b), 
even for strong interactions (see Sect. C). With many 
particles, as the strength of the interaction increases, the 
centers of mass of the clouds become increasingly coupled 
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FIG. 2: (Color online) Panel a) Time evolution of the 
spin-resolved center-of-mass XcM,a(t) of a system of twelve 
fermions in a harmonic potential. Solid lines refer to U/J — 
0.02 while dotted lines to U / J — 0. Panel b) Positions of the 
maxima of the amplitude of the center-of-mass oscillations as 
functions of time t in units of T — 2ti/uj. Different symbols 
correspond to different interaction strengths. The tiny decay 
in the non-interacting case is due to lattice effects. Solid lines 
are parabolic fits, ^CM,<x(£)| peak = X [1 - (£/tstsd) 2 ], where 
Xq is the same for all values of U/J. 



to internal degrees of freedom. If N is sufficiently large, 
5(0) becomes featureless, with a bandwidth of the order 
of uo. In this regime one has the situation of a single de- 
gree of freedom (center of mass) irreversibly transferring 
energy into a "bath" of microscopic degrees of freedom: 
accordingly, the amplitude of the oscillations decays lin- 
early in time as expected of an ordinary damped oscilla- 
tor. 

These observations imply a non-trivial crossover in the 
short-time dynamics of ^CM,cr(^)| peak as a function of 
the number of particles. In particular, as illustrated in 
Figs. 9-12 in Sec. D, we note the existence of a time scale 
t*, depending on N and U/J, below which the decay of 
the oscillation amplitudes is quadratic. The value of t* 
decreases with increasing N and increases with increas- 
ing U/J. More quantitatively, we have investigated such 
crossover by fitting numerical data at strong coupling 
with the "split-fit" formula in Eq. (Dl) of Sec. D, which 
contains tstsd and t* as fitting parameters. This equa- 



FIG. 3: (Color online) Panel a) Same as in Fig. 2a) but 
at strong coupling. Solid (dashed) lines refer to U/J = 20 
(U/J = 5) while dotted lines to U/J = oo. In the latter 
case the oscillations have twice shorter periodicity than those 
of non-interacting fermions. Note that they do not display 
an appreciable decay on the time scale of the plot. Panel b) 
Same as in Fig. 2b) but at strong coupling. Solid lines are fits 
obtained by using Eq. (Dl) in Sec. D. 



tion encodes a quadratic decay for t < t*, followed by a 
linear behavior for t > t*. From our analysis we conclude 
that the value of t* for N = 6 and 5 < U / J < 30 is much 
smaller than the period of oscillations. This explains why 
no quadratic behavior is seen in Fig. 3b). From the nu- 
merical data at strong coupling, we conclude that a linear 
decay in time of ^CM,cr(^)| peak occurs when the overlap 
between the two colliding clouds is substantial, while a 
quadratic decay takes place initially (for t < t*) when 
minor overlap occurs in the tails of the clouds. 

Our main results for the time scale tstsd associated 
with STSD are summarized in Fig. 4. Here we report the 
values of t^ sd used to produce the fits in Figs. 2b)-3b). 
We clearly see that Ts^ SD vanishes linearly in the weak- 
coupling U/J^0 limit. This has to be contrasted with 
the SD relaxation rate in a system with many degrees 
of freedom at equilibrium: the latter is quadratic in the 
coupling constant governing the strength of inter-particle 
interactions (see e.g. Ref. 7). In the strong-coupling limit 
t stsd behaves approximately like 1/U. No analytical 
results are available in this regime, even in a system with 




FIG. 4: (Color online) Empty circles represent the inverse of the short-time spin-drag time constant, Tc^ 1 ^ (in units of 1/T), as 
a function of the coupling U/J. These data have been extracted from the spin dynamics of the model (1) (quadratic fit at weak 
coupling and "split-fit" procedure at strong coupling, see Sec. D). In a wide range of coupling constants, 0.1 < U / J < 1, fittings 
like the ones in Figs. 2b)-3b) do not work. We see that t^,J sd vanishes linearly in the weak-coupling regime (long-dashed line) 
and behaves approximately like 1/U at strong coupling (solid line). The solid line is a power-law fit, i.e. 1/tstsd = A (U / J) _c \ 
with A w 1.26/T and a « 0.97. The empty triangles label t^ 1 ^ as extracted from the spin dynamics of the effective model 
(2). The short-dashed line is a power- law fit of the form 1/tstsd = B (U/J)~ p , with B w 0.7/T and ft « 0.91. While the two 
exponents a and f3 are very similar, the proportionality constants A and B are slightly different. This discrepancy may be due 
to the neglect of three-site terms [21] in Eq. (2). In the inset we show the same strong-coupling results in a log-log scale. 



many degrees of freedom at equilibrium. We emphasize 
that, for all the data at U/J ^> 1 in Fig. 4, t* is within 
the observation time of our simulations. 

To check the robustness of our conclusions at strong 
coupling, we study an effective "t- J" model [19, 20] which 
approximates (1) for U/J ^> 1: 

i ^ ' i 

(2) 

where Si = J2 a c\ a (cr a p/2)ci^ is the spin operator (a 
being a three-dimensional vector of Pauli matrices) [21]. 
When U / J ^> 1 such model is much easier to simulate 
than the original FH model. Employing Eq. (2) we have 
discovered that, for a fixed value of N, the amplitude of 
the oscillations decays quadratically in time when U/J 
is sufficiently large. This is shown in Fig. 8 in Sec. D. 
In other words, as mentioned in the Introduction, the 
quadratic dependence on time of the decay of the oscilla- 
tion amplitudes displays a reentrant behavior pertaining 
to the many-particle problem. In Fig. 4 we report the re- 
sults for the inverse STSD time constant of the model (2) 
(empty triangles), which agree qualitatively with those 
based on the full FH model. 

In summary, we studied short-time spin-density oscil- 
lations in a strongly-interacting ID few-fermion system. 
We discovered that the decay in the oscillation ampli- 
tudes goes from quadratic to linear back to quadratic in 
time as the interaction strength increases from zero to 
infinity. The inverses of the properly-defined time con- 
stants depend on the strength of inter-particle interac- 
tions in a way that was unpredictable on the basis of our 
knowledge of the same phenomenon in many-particle sys- 



tems near equilibrium. Our predictions can be tested by 
studying the damping of spin-dipole oscillations in few- 
fermion systems [15]. 
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Appendix A: Technical details on the numerical 
method 

In this Appendix we give some technical details on the 
numerical method we have used. Our simulations are 
based on the TDMRG method, which is known [13] to 
be a powerful technique for the simulation of ID systems. 

In this work we have used a matrix product state 
(MPS) representation of the wave function, enforcing sep- 
arate conservation of the number of spin-up and spin- 
down particles. The ground state at t = is found using 
the procedure described in Ref. 22 (modulo small vari- 
ations). Inversion symmetry with respect to the trap 
center, i.e. rti y(T (t) = i,o-(^) 5 is n °t enforced a priori 

but is present in the converged results (we use this fea- 
ture as one of the benchmarks of the simulations). The 
truncation step is treated in a fully "dynamical" way, 
i.e. we do not fix a maximum (or a minimum) number 
of states m for each bond link. On the other hand, we 
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FIG. 5: (Color online) Panel a) Time evolution of the spin- 
resolved center-of-mass XcM,a(t) of a system of two fermions 
in a harmonic potential. Solid lines correspond to g = 0.2, 
while dashed ones to g — 0.1. Panel b) Positions of the 
maxima of the amplitude of the center-of-mass oscillations 
as functions of time t in units of T — 2tt/oj. Different 
symbols label data corresponding to different values of the 
coupling constant g. The solid lines are fits of the form 
^CM,cr(t)| peak = X cos (a/2 £/tstsd), where X is the same 
for all values of g. 



FIG. 6: (Color online) Same as Fig. 5 but in the strong- 
coupling regime. Solid lines in panel a) correspond to g = 
150, while dashed ones correspond to g = 50. We stress that 
the fitting function used in panel b) is the same as in panel 
b) of Fig. 5, i.e. X C M,*(t)\ k = X cos (y/2 £/t S tsd). 



choose to discard states with a "small" statistical weight, 
summing up to a maximum allowed error e. This repre- 
sents the crucial parameter that controls the precision 
and the duration of our simulations. Typically, we use 
e~10 -8 — 10 -10 . We have checked that these values are 
sufficiently small by employing two different procedures: 
i) we have compared our numerical results with the ex- 
act solution that is available in the non-interacting case 
U = and ii) we have checked the accuracy in the inter- 
acting case U > by analyzing the convergence of the re- 
sults with decreasing e. With this dynamical-truncation 
procedure the maximum number of states used in the 
simulations is m ~ 1 x 10 3 — 5 x 10 3 , which is reached at 
the trap center where the time-evolution of rii y(7 is rather 
complex (as we have seen in Fig. 1). 

The time-evolution operator is treated within a Suzuki- 
Trotter expansion. Due to the presence of nearest- 
neighbor-only interactions, one can separate couplings 
on odd bonds from couplings on even bonds, thus writing 
the global Hamiltonian as the sum of two non-commuting 



terms, H = 7i e ven + ^odd- Each of the two contributions 
is the sum of commuting two-site terms. To the n-th 
order the time-evolution operator reads 

k 

e -iUAt _ e -iH even CjAt e -iH odd djAt _j_ (3(Af*+l) ^ 
J = l 

(Al) 

where the number 2k of exponentials to be multiplied 
as well as the coefficients Cj and dj depend on the or- 
der of the expansion [23]. In our simulations we em- 
ploy a sixth-order Suzuki- Trotter expansion (which in- 
terestingly enough was found to perform faster than the 
second-order one). 

Appendix B: Perturbation theory in the 
weak-coupling regime 

In the weak-coupling limit it is possible to study 
the impact of a contact repulsive interaction on the 
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FIG. 7: (Color online) Inverse of the STSD time scale tstsd 
as a function of the coupling constant g. All axes are in 
logarithmic scale. Data labeled by circles refer to the weak- 
coupling regime (lower horizontal axis), while the ones labeled 
by squares refer to the strong-coupling regime (upper horizon- 
tal axis). The solid line is a power-law fit, 1/tstsd = A g~ a 
with A w 10.4/T and a ~ 0.99. The dashed line is a power- 
law fit, 1/tstsd = C g 1 with C « 0.77/T and 7 « 0.99. 



spin-resolved center-of-mass dynamics by means of time- 
dependent perturbation theory. For the sake of sim- 
plicity, we consider a continuum model, which is more 
amenable to an analytic treatment (h = 1): 



H = J dx X>t( x )^__l 



(Bl) 



The operators i/j a (x) with a =t, I are anticommuting 
field operators for spin-up and spin-down particles, while 
Pa( x ) — ty\(%)ty(j(%) are spin-resolved density operators. 

The state of the system before the quench (t < 0) 
corresponds to N spin-up and TV spin-down particles in 
the ground states of two separate spin-resolved harmonic 
confinements. The distance 2d between the harmonic 
traps is supposed to be large enough (a few harmonic 
oscillator lengths) so that the initial state \ip(t = 0)) is 
(to a very good approximation) 



\m) 



-i p ^e iP ^ d \N f) 8> |iV 4,) , 



(B2) 



where \Na) is the non-interacting ground state of N par- 
ticles with spin a in a harmonic potential and the trans- 
lation operator P a is 

P a = -i J dx ^t(x)-^^ a (x) . (B3) 

The time evolution of the system for t > is dictated by 
the Hamiltonian (Bl). 

We now expand the two exponentials in the (small) 
parameter gt, obtaining, up to first order, the following 



XcM,a(t) = ±dcos(ut) + g 



dx 



N 



x D a (x,x',t — t')n^(x' ± 2dcos(ujt')) 



(B4) 



where the plus (minus) sign refers to a =t (cr =|), a 
denotes the spin component opposite to a, and n a (x) = 
(Na\p a (x)\Na) is the spin-resolved ground-state density 
profile. The quantity D a in the second term in the r.h.s. 
of Eq. (B4) is the density-density response function: 

D a (x,x / ,t-t / ) = -iQ(t-t') 

x (Na\[p a {x,t)rPv{x' ,t')]\Ncj) , 

(B5) 

Q(t) being the Heaviside step function. 

It can be easily shown that the density-density re- 
sponse function for the harmonic oscillator has the fol- 
lowing properties 



D a (x, x' , t-t')= D a ( — x. -x' \ t - t f ) 
D a (x, x' ', t-t')= D a (x, x' ', t — t' + T) 



(B6) 



where T — 2tt/uj. Moreover, since the first argument of 
D a (x, x' ,t — t') is integrated after being multiplied by an 
odd function [x in the first line of Eq. (B4)], only the 
antisymmetric part of the response function is needed: 



D ( ^(x,x',t-t') 



-Dl A \-x,x f ,t-t f ) 
-D^(x,-x',t-t') . (B7) 



Using the Lehmann representation one can also prove the 
following identity: 



Dl A \x,x',t) =D ( f\x,x',T/2-i) . 



(B8) 



Using these identities one can prove that the first- 
order term vanishes in correspondence of the extrema 
of XcM,a(t) (i.e., for t = utt/uj). We thus conclude that 
the envelope of XcM,a(t) decays quadratically in time, 
as shown in Fig. 2. 



Appendix C: The two-body problem 

The problem of two particles interacting with a short- 
range potential in a harmonic trap is exactly solv- 
able [24, 25]. This can indeed be reduced to a single- 
particle problem by switching to center-of-mass R = 
(xi J rX2)/2 and relative motion r = x\ — X2 coordinates. 
The first- quantized Hamiltonian in reduced units reads 

(n=i) 



Id 2 
'A OR 2 



dr 2 



l— + ^-— + y+9 5{r). (CI) 



Distances have been rescaled with the harmonic oscillator 
length aho = (mco>) -1 / 2 and energies with uj. In these 
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units g = g/ (uia^ ) is the dimensionless coupling constant 
that controls the strength of interactions. 
The initial state can be factorized as 



(x 1 ,x 2 mo)) 



(2n)i 



(C2) 



where d is the initial displacement (as in Sec. B) in units 
of aho- The center-of-mass portion of the wave func- 
tion has a trivial dynamics, while the evolution of the 
relative portion can be calculated numerically to the de- 
sired degree of accuracy by expanding the second factor 
in Eq. (C2) in the exact eigenstates [24, 25] of the Hamil- 
tonian (CI). 

In Fig. 5 we illustrate XcM,a(t) as a function of time 
t for small values of g, while in Fig. 6 we show the same 
quantity for large values of g. These plots have to be 
compared with Figs. 2-3. 

Comparing Fig. 5b) with Fig. 6b) we see that the 
time-evolution of the peak position ^CM,cr(^)| peak can 
be fitted by a functional form which is quadratic in 
time at small times both in the weak- and strong- 
coupling regimes. Indeed, to fit the data in these panels 
we have used the following function, ^CM,cr(^)| peak = 
Xq cos (y/2 £/ t stsd), which at small times reduces to 
XcM,a(t/r S TSD -> 0)| peak = X [1 - (*/t S tsd) 2 ]. This is 
at odd with the many-particle case analyzed in Figs. 2- 
3. We remind the reader that in this case the functional 
dependence of ^CM,cr(^)| peak at strong coupling on time 
is linear rather than quadratic. 

In Fig. 7 we illustrate the dependence of tstsd on the 
coupling constant g. In the weak-coupling regime we find 
1/tstsd g while in the strong-coupling regime we find 
1/tstsd oc 1/g. This is in perfect agreement with the 
numerical results shown in Fig. 4. Note also that the 
order of magnitude of tstsd/^ in Fig. 7 is the same as 
in Fig. 4. 



Appendix D: Split-fit procedure for the 
linear-to-quadratic crossover 

The mapping of the Fermi-Hubbard model into the 
effective "t-J" model (Eq. (2)) allows to explore a wider 
range of values of the coupling constant U/J. Typical 
results for the spin dynamics of this model are shown 
in Fig. 8. Notice that for very large couplings (U/J ~ 
100 — 200) the data for the peak position ^CM,cr(^)| peak 
are not well fitted by a linear function of time, and a 
quadratic term becomes non-neglibible. Therefore, the 
peculiar linear decay of the oscillation amplitude that is 
visible in Fig. 3b) seems to disappear (or, better, to take 
place at later times) for very strong couplings. 

We have thus carried out a detailed numerical analysis 
of the linear-to- quadratic crossover, both as a function 
of the number of fermions with a given spin N and of 
the coupling constant U/J, for both Fermi-Hubbard and 



a) 



b 



b) 




FIG. 8: (Color online) Same as in Figs. 2-3 but for the effec- 
tive "£-J" model defined in Eq. (2). In panel a) the solid 
lines correspond to U/J = 20, while the dotted ones to 
U/J = oo. In panel b) the solid lines are quadratic fits of 
the form XcM,a \ (t) pea k = Xo(l + bt + at 2 ). We indicated 
with red curves (data for JJ / J — 100, 200) the cases in which 
the quadratic term in the fit is not negligible. 



"t-J" models. Our procedure is based on the key ob- 
servation that the quadratic behavior takes place from 
t = to a crossover time t = t*(N,U/J), after which 
the oscillations decay linearly. This is clearly seen in 
Fig. 9, where we illustrate the time evolution of the spin- 
resolved center-of-mass ^CM,cr(^)| peak f° r different values 
of N and U/J and for the effective "t-J" model. 

As a matter of fact, we have adopted a systematic 
"split-fit" procedure by fitting our data at strong cou- 
pling with the following function: 



^CM,aW| peak 



Xoll 



Xn 



t A 



2t*TSTSD / ' 
t t* 



t < t* 



TSTSD 2TSTSD 



t > t* 

(Dl) 

which yields both tstsd, the short-time spin-drag time 
constant defined in Sec. Ill, and t*(N,U / J). The r.h.s. 
of Eq. (Dl) has been expressly written in such a way to 
guarantee that ^CM,cr(^)| peak is a continuous function of 
time t, with continuous first derivative, at t = t*. 

From Fig. 9 we clearly see that, as N is increased 
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or U/J is decreased, the time scale t*, below which 
^CM,cr(£)| peak is well fitted by a parabolic function, de- 
creases. The values of t* and tstsd extracted from this 
procedure are reported in Fig. 10 and also in Fig. 4. 



We have repeated the same numerical analysis also for 
the Fermi-Hubbard model. The results are reported in 
Figs. 11-12 and also in Figs. 3-4. 
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FIG. 9: (Color online) Time evolution of the spin-resolved center-of-mass XcM,a(t)\ peak of a system of 2N fermions in a 
harmonic potential. From top to bottom N decreases from 7 to 1. The three different panels refer to different values of the 
coupling constant U/J. The fitting function is given in Eq. (Dl). The initial quadratic decay for < t < t* is shown by 
solid red lines, while the subsequent linear decay for t > t* by dashed blue lines. Data in this figure have been produced by 
employing the effective "t-J" model - Eq. (2). 
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FIG. 10: (Color online) Dependence of t* [panel a)] and tstsd [panel b)] on U/J and N for the effective "t-J" model. These 
data have been extracted by applying the split-fit procedure (Dl) to the data reported in Fig. 9. 
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FIG. 11: (Color online) Same as in Fig. 9, but only for U/J — 10. Panel a) refers to the Fermi-Hubbard model, while panel 
b) refers to the effective "£-J" model. Note the different horizontal scale: the simulation of long-time spin dynamics of the 
Fermi-Hubbard model is more difficult. 




FIG. 12: (Color online) Same as in Fig. 10, but only for U/J — 10. Here we compare results for the Fermi-Hubbard model 
[Fig. 11a)] with those for the effective "t-J" model [Fig. lib)]. Note that the values of t* for the two models are practically 
identical. 



